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ABSTRACT 

Motivated by recent progress in the statistical modeling of quasar variability, we develop a new ap- 
proach to measuring emission-line reverberation lags to estimate the size of broad-line regions (BLRs) 
in active galactic nuclei. Assuming that all emission-line light curves are scaled, smoothed, and 
displaced versions of the continuum, this alternative approach fits the light curves directly using a 
damped random walk model and aligns them to recover the time lag and its statistical confidence 
limits. We introduce the mathematical formalism of this approach and demonstrate its ability to cope 
with some of the problems for traditional methods, such as irregular sampling, correlated errors, and 
seasonal gaps. We redetermine the lags for 87 emission lines in 31 quasars and reassess the BLR 
size-luminosity relationship using 60 H/3 lags. We confirm the general results from the traditional 
cross-correlation methods, with a few exceptions. Our method, however, also supports a broad range 
of extensions. In particular, it can simultaneously fit multiple lines and continuum light curves which 
improves the lag estimate for the lines and provides estimates of the error correlations between them. 
Determining these correlations is of particular importance for interpreting emission- line velocity-delay 
maps. We can also include parameters for luminosity-dependent lags or line responses. We use this 
to detect the scaling of the BLR size with continuum luminosity in NGC 5548. 
Subject headings: galaxies: active — galaxies: nuclei — galaxies: Seyfert — quasars: general 



1. INTRODUCTION 

While it is widely accepted that the enormous lu- 
minosities of active galactic nuclei (AGNs) are at- 
tributable to accretion of matter onto supermassive black 
holes (BH), detailed studies are extremely challenging 
on account of the small angular scales of the regions in- 
volved in the accretion process. Direct probes of the 
sub-microarcsecond structure of AGNs has been there- 
fore limited to VLBI studies of the radio-emitting re- 
gions, gravitationa l microlensing stud ies of the accretion 
disk (see review by [Wambsganss 2006) and reverberation 
mapping of the broa d- line regions dBlandford fc McKed 
Il982t IPetersoiil 11993ft . The technique of reverberation 
mapping (a.k.a. echo mapping) exploits the light travel 
time between the central engine and the broad-line re- 
gion (BLR) to deduce the structure of the BLR (sec 
iPetersonl 120011 for a tutorial). The continuum radiation 
from the accretion disk photoionizes gas clouds near the 
AGN to produce broad emission lines, thus encoding the 
geometry and kinematics of the clouds (jOster brocldl 1 9891 : 
lPetersorJ[l997l lKrorildll999t ). The physical ansatz for re- 
verberation mapping is straightforward: 

1. The continuum emission of the quasar shows 
(stochastic) variability that drives emission-line 
variations after a light travel-time delay. 

2. The unobservable ionizing UV continuum that 
drives the emission lines is simply related to the 
observable satellite UV or optical continuum (i.e., 
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the pattern and phase of variations are closely cor- 
related) . 

3. The light-travel time is the most important time 
scale; specifically, the local emission-line response 
time to continuum changes is assumed to be instan- 
taneous and the dynamical time scale of the BLR 
is much larger than the light-travel time across it. 

The relationship between the observables, continuum 
light curve s c (t) and the emission- line light curve si(t, V) 
where V is the line-of-sight velocity, is taken to be 

ai(t,V) = J dr^{r,V)s c {t-r), (1) 

where \P(r, V) is known as the "transfer function" or 
"velocity-delay map." In reality, the relationship be- 
tween the continuum and emission-line variations can be 
non-linear, but the amplitude of variation on reverbera- 
tion time scales is sufficiently small that the linear ap- 
proximation seems to be justified. Inspection of Equation 
([1]) shows that ^(r, V) is the observed response of the 
broad emission-line region to a delta-function continuum 
outburst, mapped into the observable quantities time de- 
lay r and line-of-sight velocity V. The data requirements 
for successfu l recovery of the tr ansfer function are quite 
demanding ([Home et al.l 120041) and consequently most 
efforts to date have concentrated on measuring only the 
total emission-line response to continuum variations. The 
transfer equation (cq. [T]) then becomes 

si(t) = J dr*(r)s c (t-r), (2) 

where ^(t) = J ^(r, V) dV is variously known as the 
"one-dimensional transfer function" (so that ty(t,V) is 
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the "two-dimensional transfer function" or the "delay 
map"). For the remainder of this paper, we will refer to 
^f(t) simply as the "transfer function." In most investi- 
gations to date, it is the mean response time or "lag" (r) 
that one tries to measure, generally by cross-correlation 
of the continuum and emission-line light curves, as we 
discuss further below. The importance of measuring the 
emission-line lag is two-fold: first, (r) yields a charac- 
teristic physical scale for emission of a particular line, 
R = c(r), and this can be combined with some measure 
of the emission-line Doppler width AV to obtain an es- 
timate of the central black hole mass. Assuming that 
gravity is the dominant force on the line-emitting gas, 
the virial equation for the central black hole mass is 



Mi 



Bl] 



fAV 2 R 
G ' 



(3) 



where G is the gravitational constant and / is a dimen- 
sionless factor of order unity that depends on the geome- 
try, velocity field, and inclination of the BLR. We note in 
passing that there is currently an active debate about the 
relative importance of radiation pressure on the BLR gas 
and ho w this affects reverberation-based mass measure- 
mcnts dMarconi et alj|2008t iNetzerl 120091 i Marconi et all 
120091: iNetzer fc Marzianil 120101 ) . While the possible role 
of radiation pressure in measurement of black hole masses 
is an important issue, it has no direct bearing on the 
present discussion, which is about measuring time delays. 
Similarly, there is still active discussi on about the mean 
value of the scaling factor (/) (e.g.. lOnken et al. 1200 
lLabita et al.1[200l I Woo et al.ll20igflGraham et al.ll20TTl 



that is beyond the scope of this contribution. Second, 
reverberation studies have established a tight empirical 
relationship betwee n the BLR size and the AGN con- 
tinuu m luminosity dKaspi et alJll996t iKaspi et alJl2000L 
120051 iBentz et alj 120061 12009h that allows us to use the 
luminosity as a surrogate for the BLR size in eq. ^ and 
thus estimate the masses of black holes in AGNs from in- 



dividual spectra (IWandel et al.l 


19991: iVestertraardl 120021 


McLure & Jarvis! 20021 I2004 


Vestereaard & Peterson! 


20061 IKollmeier et al.l 120061: IS] 


ren et al.l 120081). This 



allows us to explo re BH properties and evolution 
with redshift (e.g.. IKollmeier et all [20061 iPeng et all 
20061: iHopkins fe Hernquistl 120061: iShankar et al l 120091: 



Steinhardt fe Elvisl 120101 i IKelly et al.l 120101 ) . thus pro- 
viding valuable insights into the mystery of black hole 
growth and its connection to galaxy evolution at high 
redshift, where the quasar population is evolving dra- 
matically. The potential for obtaining simple esti- 
mates for the masses of black holes in quasars pro- 
vide a means of exploring the correlations between the 
BH mass and global properties of their host galax- 
ies such as the bulge luminosity ( Mrm-Lhu]^ rela- 
tionshi p: iKormendv fe Richstonelll995l ; iMagorrian et al.l 
119981: IBentz et all |2009|) and bu lge stellar velocity 
dispe rsion (Mbh~0* relationship : iFerrarese fc MerrittJ 



20001 Gebhardt et al . 2000a.b: iTremaine et al.l 120021: 



Ferrarese et al.1 120011: lOnken et alJ 120041 : iNelson et ail 



20041: iGultekin et alJ 120091 ) both locally and potentially 



over cosmic time. 

As a practical problem in aperiodic time-series data 
analysis, reverberation mapping requires high-fidelity 
spectroscopic monitoring of the continuum and emission- 



line variations fo r a duration long compared to the 
emission-line lag dHorne et al.l 120041) . which is observed 
to range from hours to a year or more, depending on 
the luminosity of the AGN and the cosmic time dila- 
tion at its redshift. Emission-line lags have been mea- 
sured for more than three dozen AGNs by cross corre- 
lation of the continuum and emission-line light curves. 
The particular challenge of dealing with reverberation 
time series is that they are generally irregularly sam- 
pled for various reasons, including unfavorable weather 
and, for higher-luminosity objects with larger lags, an- 
nual conjunctions with Sun that cause seasonal gaps in 
the observed time series. In practice, two methodolo- 
gies have been widely employed to deal with unevenly 
sampled data. The first method is to interpolate be- 
tween real data points to obtain a regular sampling 
grid for computation of the cross-co rrelation function 
(CCF ) as a function of time delay r flGaskell fc Sparkd 
1981 IGaskell fc PetersoriH987t IWhite fc Petersonlll99T 



Peterson et al. 1 119981: IWelshl 11999b IPeterson" et al. 

The second me thod, the discrete corre lation function 
(DCF) method (|Edelson fc Krolild fl98l . bins the data 
over discrete time intervals on which the data are rea- 
sonably well-sampled and a correlation coefficient is 
computed for the time delays between each pair of 
continuum/cmission-line time bins. A var iant on this is 
the ^-transformed DCF (I Alcxand"erl 119971 ) which varies 
the width of the time bins to better distribute the data 
points among the time bins. IWhite fc Peterso"nl ()1994f ) 
show that when common assumptions and normaliza- 
tions are used, the interpolation CCF method and the 
DCF method give similar results. However, as the time- 
sampling becomes sparser, the interpolation method sig- 
nificantly outperforms the DCF method as long as inter- 
polation of the light curves (usually linear in practice) 
remains a reasonable assumption^]. 

Presumably even more accurate lags could be mea- 
sured given more realistic modeling of the continuum 
behavior between real measurements of the continuum 
and emission-line fluxes. This now seems to b e a rea l 
possibility given the recent work of IKelly et all (|2009l) . 
who find that quasar variability can be well described 
by a damped random walk. By applying the variability 
model to the light curves of known q uasars and com- 
paring them to other variable sources, iKozlowski et all 
(|2010l ) show that quasars occupy a very distinctive re- 
gion in the model parameter space of tim e scale and vari- 
ability amplitude. iMacLeod et al.l (|2010l ) then apply the 
model to ~ 9,000 spectroscopically identified quasars in 
SDSS Stripe 82. They confirm that the model can de- 

3 Consider, as an example, the UV and optical monitioring 
campaign on NG C 5548 undertaken with the International Ul- 
traviolet Explor er (Clavcl et al. 1991) and ground-based telescopes 
(IPeterson et al.H199H) in 1989. The UV data were sampled at ap- 
proximately regular 4- day intervals. Analysis of these data using 
the interpolation CCF (Peterson & Wandcl 1999) revealed for the 
first time a "virial relationship" between the emission-line lags and 
line widths (i.e., (t) oc AV~ 2 ), thus providing an empirical jus- 
tification for using Equation 10 to estimate the blac k hole mass. 
Anal ysis of these same data with the DCF method (Krolik et al. 
[1991J) obscured this result at least in part because of "discretization 
noise" introduced by the DCF: because of the regular 4-day sam- 
pling, the smallest usable time bin for the DCF was also four days, 
resulting in lag measurements that were integer multiples of 4 days, 
significantly reducing the time resolution of the lag measurements 
and smearing out the (t)-AV anticorrelation. 



Measuring Reverberation Lags 



3 



scribe quasar variability well and they explore the corre- 
lation of the variability parameters with other properties 
of quasars such as wavelength, luminosity, BH mass, and 
Eddington ratios in detail. Importantly for reverbera- 
tion studies, the formalism is able to statistically predict 
the value of light curve at an unmeasured time based 
on the overall statistical properties of the light curve. It 
provides a well-defined statistical model for interpolating 
light curves and can do appropriate statistical averages 
over the uncertainties in the model predictions. 

Given a complete statistical framework for describing 
the continuum variability, and the overall ansatz that 
emission-line variability is a scaled and smoothed version 
of the continuum, we can build an alternative approach 
to measuring reverb eration lags, aspec t s of w hich were 
previously noted by iRvbicki fc Klevnal (|1994[ ) . Among 
the advantages of this approach are: 

1. It not only interpolates between data points, but 
also self-consistently estimates and includes the un- 
certainties in the interpolation. 

2. It can separate light curve means, trends, and sys- 
tematic errors in flux calibration from variability 
signals and mcansurcmcnt noise in a self-consistent 
way. 

3. Correlated errors can be treated naturally. 

4. Lags of multiple emission lines and their covari- 
ances can be derived simultaneously. 

5. It provides statistical confidence limits on the lag 
estimates as well as other parameters. 

We describe the methodology of our approach in detail 
in fj2j In Sj3l we present the statistical process model for 
the continuum light curves. We briefly describe our data 
set and apply this method to the estimate of H/3 lags 
in iJU We further show how the method can address the 
problem of correlated errors in $5] and how it can be used 
to improve lag estimates, particularly in the presence of 
seasonal gaps, by fitting multiple lines simultaneously 
in £j6j We also fit the -Rblr - relationship using H/3 
lags determined by our method in $7] In EjSj we add a 
luminosity dependence to the lag and solve for the lag- 
luminosity relationship of NGC 5548. We summarize 
our main findings and discuss future applications and 
expansions of our approach in fj9] 

2. METHODOLOGY 



iPress et all (fl99l and IRvbicki fc Presi (|l99l devel- 
oped a method to sta tistically analyze irregula rly sam- 
pled light curves, and IRvbicki fc Klevnal (|1994h applied 
the variant we now consider to four seasons of optical re- 
verberation data on NGC 5548. Here we reintroduce this 
approach, which we have named "Stochastic Process Es- 
timation for AGN Reverberation (SPEAB0)," with sev- 
eral modest changes in algorithm and a broad range of 
new applications. 

Except for the transfer fu nction ^>(t), our notatio n 
is chosen for comparison with IRvbicki fc Klevnal (|1994h . 

4 http: //www. astronomy . Ohio- state . edu/~yingzu/spear .html 



We start with a model process driving the continuum 
s c (t) that has a covariance between times ti and tj of 



{s c (ti)s c {tj)) = a 2 exp(-|ij - ^ |/r d ). 



(4) 



We adopt here an exponential co variance ma t rix fo r 
concreteness, since we kno w from iKellv et al.l (12009ft . 
iKozlowski et all (|2010j) and iMacLeod et al.l (|2010t) that 
quasar light curves are well modeled by this process. 
Physically, the model corresponds to a random walk de- 
scribed by an amplitude a 2 = <7 2 Td/2 on long time scales 
and an exponential damping time scale th, where a and 
Td are used as our model parameters. IRvbicki fc Press! 
(|1992f) estimated the covariance matrix based on the 
structure function of the continuum light curve, while 
here we adopt a specific parametrized model that will be 
optimized as part of the analysis. 

Slightly rewriting Equat ion ([2t for conveni e nce a nd to 
facilitate comparison with IRvbicki fc Klevnal (|1994l ) , the 
light curve of a line is 



si(t) = / dt'^(t-t')s c (t). 



(5) 



Since the lines and continuum are related by the transfer 
function, we can also determine the covariance between 
the line and continuum 

(si(U)s c (tj)) = [ dt'y(U-t')( Sc (t')sc(t 3 )), (6) 



between the line and itself 



(siiUUit,)) = / dt'dt"*(t i -t , )*(t j -t")(s c (t')s c (t")), 



and between two different lines 



(7) 



( S i(tM(tj)) = J dt'dt"*(u-t')m'(t 3 -t")( Sc (t')sc{t")). 

If the light curve of the line is divided into velocity 
bins 5V , then there is a transfer function for each bin 
\&(t — t' , V) and we can compute all the expected covari- 
ances between the light curves. For convenience, let s 
be a vector comprised of all the light curves, both line 
and continuum, and S — (ss) be the covariance matrix 
between all the elements of s. By definition, in Gaussian 
statistics the probability of the light curve is simply 



P(s) oc |sr i/2 



exp 



s T S- 



(9) 



We do not measure the actual light curve, but some re- 
alization of it, y = s + n + Lq, in which there is mea- 
surement error n, whose probability distribution is 



P(n) oc \N\ 1/2 exp ( - 



(10) 



where N = (nn) is the covariance matrix of the noise. 
Note that nothing requires N to be diagonal, so there is 
no formal difficulty to including covariances in the noise 
between the line and continuum. 

In defining y, we have also allowed for the simultaneous 
fitting of a general trend defined by a response matrix L 
and a set of linear coefficients q. In particular, we use this 
to fit and remove separate means from the light curves. 
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X 2 /dof 

Fig. 1. — Distribution of x 2 per degree of freedom for the con- 
tinuum fits. The solid histogram is the x 2 /dof distribution of our 
stochastic model, while the dashed one shows the distribution ex- 
pected for models with correctly estimated Gaussian uncertainties. 
The dotted histogram is the x 2 /dof distribution of the joint model 
of the continuum and H/3 light curves from [|4] 

In this application to a model with two light curves, L is 
a 2 x A" matrix with entries of (1,0) for the continuum 
data points and (0, 1) for the line data points, where K 
is the total number of data points. The linear parame- 
ters are a very general tool. For example, separate linear 
trends would be removed with a 4 x K matrix with en- 
tries of (1, ti, 0, 0) for continuum epoch ti and (0, 0, 1, tj) 
for line epoch tj. Two sources of data with potentially 
different but constant levels of contamination from the 
host galaxy can be reconciled by using different means 
for each line and continuum data source, corresponding 
to a 4 x K matrix with entries of (1, 0, 0, 0) for the first 
continuum source, (0,1,0,0) for the second continuum 
source, (0,0, 1,0) for the first line source and (0,0,0, 1) 
for the second line source. Unlike current approaches fo- 
cused on cross-correlation functions, the uncertainties in 
these linear parameters arc fully incorporated into the 
uncertainties in any other parameter estimate. 

Given these definitions, the probability of the data y 
given the linear coefficients q, the intrinsic light curves 
s, and any other parameters of the model p is 



P(y|q,s,p) cx |£W| 



-1/2 



d n n d n s S (y - (s + n + Lq)) exp 



B T S- 1 B + n T N- 1 n 



(11) 



After evaluating the Dirac delta function, we "complete 
the squares" in the exponential with respect to both the 
unknown intrinsic source variability s and the linear co- 
efficients q. This exercise determines our best estimate 
for the intrinsic variability 

§ = SC- 1 (y-£q) (12) 

and the linear coefficients 

q = (L T C- 1 L)- 1 i T C- 1 y = C^C^y (13) 
where C = S + N is the overall covariance matrix of the 



data and C q = (L T C X L) l . With these definitions we 
can factor the argument of the exponential into 

P(y|q,s,p) k\SN\- 1/2 

As T (S- 1 +N~ 1 )As _ Aq T C- 1 Aq _ y T C^y 
2 2 2 

(14) ' 



exp 



where 



CJ 1 = C- 1 - C^LCqL' C 



T/~i-1 



(15) 



is the component of C that is orthogonal to the fitted 
linear functions, the variances in the linear parameters 
are 

(16) 



(Atf) = (L T C- 1 L)-^C, 



As = s — § and Aq = q q. We can marginalize the 
probability over the light curve s and the linear param- 
eters q under the assumption of uniform priors for these 
variables to find that 



P(y|p) k£. 

= \S + N\- 1/2 \L T C- 1 L\ 



-1/2 



exp 



(17) 



where C represents the likelihood function we are to max- 
imize, and the remaining parameters p are those de- 
scribing the process (Equation [4} and the transfer func- 
tions. The term in the exponent, y T Cj y, is the gener- 
alized x 2 that we present throughout the paper. While 
this treatment o f line a r parameters was includ ed by 
iRvbicki fc Press! (|l99l . iRvbicki fc Klevnal (|l99l chose 
to subtract fixed means rather than marginalizing over 
them as part of the analysis as we do here. The variance 
in the estimate for the mean light curve is 



(As^ 



S — S^CiS. 



(18) 



We can estimate the light curve s(t) at any unmeasured 
time using the same formalism. The simplest means of 
doing so is simply to pad the data vector y^ with ad- 
ditional fake points y/ that have infinite measurement 
uncertainties in the sense that A -1 — > for these points. 
After appropriately partitioning the matrices, the esti- 
mate of the light curve at the unmeasured points is 

s f = Sfd(Sdd + N dd )- 1 y d (19) 

with variance relative to the true light curve of 



{As}) = S ff - S fd (S dd + N^Sdf. 



(20) 



where S dd , Sff, Sf d and S d f are the data-data, fake- 
fake, fake-data and data-fake covariance matrices of the 
process and N dd is the noise matrix of the data. The 
inclusion of the fake points has no effect on the expected 
results for the measured data points. 

Just to re-emp hasize the point, t hi s form alism was 
first outlin e d by IRvbicki fc Klevnal ( 1994fl based on 
IPress et al.1 (|1992D and lRvbicki fc Press! (| 199% . We have 
refined it slightly to use a specific process model, to op- 
timize the parameters of that model and to include the 
means of the light curves as parameters that are auto- 
matically marginalized. Unfortunately, we will not be 
able to use the fast implementation of this method for 
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Fig. 2. — Continuum mo dels . The solid line shows the expected mean source light curve s ( Equation I12[ l and the dashed line shows the 
expected spread (Equation 1181) of light curves about the mean consistent with the data. An individual light curve realization consistent 
with the data (see Equation l21H will show more structure than this mean light curve and have excursions outside the dashed line consistent 
with the estimated variance. 



expon ential covariance matrices from iRvbicki fc Press! 
(|1995f) . because the inclusion of the transfer functions 
means that S is not a simple exponential covariance ma- 
trix and hence does not have a simple, tridiagonal inverse 
for the fast method. 

We can, however, use the fast methods for generating 
simulated light curves. In particular, we are interested in 
light curves constrained t o resemble the co n tinuum light 
curve. As discussed by IRvbicki &: Press! (|1992f ) . such 
a light curve is simply the estimated mean light curve 
given by Equation (|12j) with an added random compo- 
nent that has the covari ance matrix Q = (S^ 1 +-/V~ 1 )~ 1 . 
IRvbicki fc Press! ()1992l ) suggest deter mining the eigen- 
modes of Q which are then the independent "normal" 
modes that can be added to the mean light curve to pro- 
duce a random realization constrained by the continuum 
light curve. This is computationally expensive. Instead, 
we note that if we Cholesky decompose Q = M T M, 



where M is an upper triangular matrix, and define the 
random component of the light curve by u — Mv where 
r is a vector of zero-mean, unit-dispersion Gaussian ran- 
dom deviates, that 



(uu T ) = M(yy t )M t = MM T = Q T = Q 



(21) 



since the covariance matrix (rr T ) of the Gaussian devi- 
ates is simply the identity matrix and Q is symmetric. 
Since Q _1 is a tridiagonal matrix given the exponential 
covariance matrix and a diagonal noise matrix, we can 
generate very high dimension u that can be convolved 
with the transfer function to produce a simulated line 
light curve in O(K) operations rather than the 0(K 3 ) 
needed following the eigenmode approach. 

The original application of the method by I Press et al.l 
(|1992t) was to cross correlate the light curves of two im- 
ages of a lensed quasar in order to estimate the time 
delay between them. While this was not discussed in 
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terms of transfer functions, it does correspond to a trans- 
fer function of the form ^(ti — tj) = S(U — (tj + At)), 
maki ng the second ligh t curve a lagged version of the 
first. iPress et al.1 (|1992j ) also treated the parameters cor- 
responding the process as fixed parameters, derived by 
fitting a power law to the structure functions of the light 
curve. It is likely that some combination of neglecting 
uncertainties in the process model or covariance s in the 
errors of the light curves led IPress et all (|1992f) to ob- 
tain an incorrect estimate of the time delay despite the 
elegance of the approach. 

In Rybicki & Kleyna's (1994) expansion of the method 
to reverberation mapping, they used rising and falling 
sawtooth and isosceles triangle transfer functions, finding 
little difference between the results or ability to discrimi- 
nate between them. Thus, for this initial reconnaissance, 
we will simply use a top-hat (rectangular function) for 
the transfer function, 

$(t-t') = A(t 2 -ti)~ l for ti<t-t'<t 2 (22) 

which has a mean lag of (t) = (t\ + t 2 )/2 and a temporal 
width of At = t 2 — t\. The necessary integrals for Equa- 
tions ©, (JTJ), and §E§ are all analytic (see the Appendix) 
and the model includes the limits of a delta function as 
At — » and a uniform thin shell as t\ — > 0. The scal- 
ing coefficient A determines the line response for a given 
change in the continuum (i.e., the responsivity of BLR 
clouds), but for present purposes we will largely view it 
as a nuisance variable. 

We use the amoeba minimization method (| Press et al.l 
Il992f) to optimize the solution a nd then either a Monte 
Carlo Markov Chain (MCMC, iMetropolis et al.l 119531: 
iHastingsl Il970) or optimization over a grid to estimate 
parameter uncertainties. We carry out the analysis in 
two phases. We first analyze the continuum light curve 
on its own, using logarithmic priors for Td and a to deter- 
mine the range of the variability process parameters con- 
sistent with the continuum light curve. The logarithmic 
prior on Td essentially penalizes values that deviate from 
the median sampling intervals to avoid both unphysi- 
cally large Td and a second class of solutions of Td — > 0, 
when all data are completely uncorrelated and the model 
simply uses a to broaden the uncertainties until obtain- 
ing an acceptable fit. Then we do the joint analysis of 
the continuum and the lines using Gaussian priors for Td 
and <t determined from the analysis of the continuum in 
isolation. In detail, we take the results of the MCMC 
analysis of the continuum and used uncorrelated priors 
on In Td and In a (which is conservative) , where the prior 
for each variable was centered at the median value with 
the Gaussian width chosen to match the upper and lower 
la confidence regions. We then used uniform priors for 
A, t\ and t 2 . 

The reason for using the continuum to define a stronger 
Gaussian prior on the process variable before carrying 
out the joint analysis is to eliminate the aforementioned 
second class of solutions of Td — > that could poten- 
tially bias our lag estimates. This secondary solution 
always exists at some level because of the finite temporal 
sampling. For modeling the continuums, we are only an- 
alyzing cases with significant variability, so this is not an 
issue for the individual light curves. However, in the joint 
analysis, if we fit the line and continuum light curves si- 
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Fig. 3. — Rest-frame damping timescale Td of the continuum light 
curves as a function of optical luminosity. The uncertainties in Td 
are the ±lcr range. 

multaneously at the wrong lag, the optimal solution will 
be to let Td — > since there are then no correlations 
between data points. Physically, it made more sense to 
consider only the ranges for the process variables Td and 
er that were statistically consistent with the continuum 
variability. 

3. THE STATISTICAL PROCESS MODEL OF THE 
CONTINUUM 

This approach depends on using a statistical model for 
the variability process of the continuum in order to opti- 
mally model the underlying light curve of the continuum. 
He re we use the expon ential covariance matrix suggested 
by IKellv et al.l (120091) . alt hough it was also introduced 
bv lRvbicki fc Press! (|1995f ) to enable a fast version of the 
SPEAR approach. Physically, the exponential covariance 
matrix in Equation |4] corresponds to a damped random 
walk with an amplitude scale a and a damping time scale 
Td . On long time scales the variance of the light curve is 
o- fTd/2) 1 / 2 an d on short time scales it is ay/i. 

Kelly et all (|2009h use this to model the light curves 
of 100 quasars, including some of the objects we will 
consider here, using a light curve fore casting approach 
to est imate the proc ess paramet e rs. | Kozlowski et al.l 
(|2010ft show how the IKellv et al] (<2009D approach can 
be derived from the SPEAR approach and demonstrated 
that forecasting is less statistically optimal for parame- 
ter estimation than using the complete light curve mod- 
eling method of SPEAR, and then applied the pro- 
cess inodel__aiidthe__SPEAR method to the OGLE- 
III (lUdalski et all l2008h light curves of ~ 2500 mid- 
infrared-selected quasars behin d the Magellanic Clouds 
(jKozlowski fc Kochanekl I2009D . They confirm that 
the damped random-walk model describes quasar light 
curves well, and that quasars occupy a well-defined re- 
gio n of Td-g parameter sp ace. This is further confirmed 
by iMacLeod et al.l (|2010f ) . who used this approach to 
model 9,000 SDSS quasars to examine the correlations 
of a and Td with other quasar properties. 

Unlike the previous papers, we fit flux rather than mag- 
nitude light curves because the line flux is more closely 
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Fig. 4. — Comparison of rest-frame H/3 time lags from the CCF and the SPEAR methods. Green triangles, red squares and blue circles 
were used for the PG, NGC and other objects, respectively, with ±lcr error bars indicated on both estimates. The labeled points linked by 
dashed vertical lines are objects having multiple lag solutions and the filled symbol is the higher likelihood solution. The two intersecting 
stripes indicate the region where the solutions from both methods may be false due to the seasonal gap (140-200 days, with time dilation). 



related to the continuum flux than to the continuum 
magnitude. Thus, we start by examining how well the 
damped random-walk process models the 60 continuum 
flux light curves for the 31 systems we consider in Sj4j 
Figure [T] shows the distribution of the x 2 P cr degree 
of freedom for the best-fit models of all the continuum 
light curves we consider. Since half of the continuum 
light curves in our sample have less than 50 data points, 
the expected x 2 /dof distribution is broader than that 
of the OGLE light curves (~ 500 points) considered by 
iKozlowski et all (|2010f ). Nevertheless, the x 2 /dof distri- 
bution indicates that the statistical process model pro- 
vides a reasonable fit to the light curves. The fact that 
the distribution is narrower than expected for correctly 
estimated Gaussian uncertainties suggests that the re- 
ported photometric errors are somewhat larger than the 
true uncertainties, or that there has been some pruning 
of outliers from the light curves. 



Figure [2] shows three examples of modeled continuum 
light curves interpolated and extrapolated from Equa- 
tion (fl"2"|) and their uncertainties from Equation (fT5)l , as 
well as the observed light curve. The estimated light 
curve at time t is in essence a weighted average over data 
points within the damping time \t — t'\ < Td that bal- 
ances the variance expected on those time scales due to 
the process against the uncertainties in the data point to 
determine how closely the model light curve approaches 
a particular data point. Far from any data points, the 
model returns to the light curve mean on the time scale 
Td. Remember, however, that Equation (|12j) is an esti- 
mate for the average of all possible light curves that could 
be drawn from the process that would be consistent with 
the data — a particular realization of s uch a light curve 
would show additional structure (see iRvbicki fc Press! 
EUg) . The "error snake" surrounding the model light 
curve is the variance in these possible light curves. Near 
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Fig. 5. — Sensitivity of the lag estimates to the noise correlation coefficient r between the H/3 and the continuum light curves of PG 0844. 
The left top panel shows the dependence of the lag on the correlation coefficient r. The left bottom panel shows the corresponding change 
in the likelihood function with r at the best-fit lag. In these panels, the blue triangle, green circle, and red square mark the results for 
r = —1,0, +1, respectively, and the dotted line indicates the 3cr limit of the likelihood function. The right panel compares the lag likelihood 
distribution for these 3 cases: r = —1 (blue dotted curve), (green dashed), and +1 (red solid), respectively. The dashed lines in the two 
right panels indicate the position of the best-fit lag, which is almost the same for all 3 cases. 



data points, its width approaches that of the measure- 
ment errors and then grows as the distance At from any 
data point increases. The variance from the process ini- 
tially increases as erIAtI 1 / 2 , but then saturates at the 
overall process variance once |A£| 3> Td. Thus, in the 
extrapolated regions we see the model light curve be- 
comes a constant and the error snake expands and then 
becomes constant. 

The three objects shown in Figure [2] represent three 
typical levels of light-curve sampling quality for the ob- 
jects we consider. Generally, the light c urves of the 
Palom ar-Green (PG) quasars obtained by IKaspi et al.l 
(|2000t) were sampled every 1-4 months over a baseline as 
long as 7.5 yr, as opposed to most of the low-luminosity 
Seyfert 1 AGNs that were more densely sampled over 
shorter baselines. The rest of the sam ple mainly con- 
sists of nearby bright Seyfert galaxies (jPeterson et al.l 
119981) whose lig ht curves are s parsel y sampled over a 
short baseline. IPeterson et al.l (|2004[) discuss the data 
in detail. In addition, we also include new light curves 
from a recent high sampling rate, multi-month rever- 
ber ation mapping camp aign on six local Seyfert galax- 
ies (jDennev et al.ll2010l) . 

We expect the damping timescales Td to show cor- 
relations with the physical characteristics of the accre- 
tion disk such as the ma ss of the centr al black hole, 
and t he AGN luminosity (|Petersonl [20081) . Kelly et al.l 
(|2009f) demonstrated this scaling relationship between 
Td and Lagn b y performing a lin e ar reg ression of Td 
on £agn, while ICollier fc Peterson! ((20011 ) also found a 
positive correlation between the characteristic timescale 
and black hole masses. Their characteristic timescale, 
which is defined by the timescale where the structure 
functions flattened, is roughly equivalent to Td. Fig- 
ure [3] shows that the more luminous central engines have 
longer exp onential damping t imescales, as we would ex- 
pect from IKellv et ail (|2009D . up to any minor differ- 



ences fronT_Jh^tin£fluxes_ra,ther than magnitudes. Note 
that iMacLeod et ail (|2010l) argue that the dependence 
of Td on black hole mass Mbh is the real driver of the 
correlation between Td and luminosity. We can use these 
correlations to estimate Td for sources lacking sufficiently 
good light curves. 

4. ESTIMATING EMISSION-LINE LAGS 

As our first application of the SPEAR method we 
recompute the lags of 101 emission-line light curves 
for 31 objects in t he literature (the compilation of 
Peterson et al.l I2004 | with the add ition of data from 
Bentz et al l 120061 iGrier et al.l 120081 and iDennev et al.l 
2006 LI2010D . We carried out the analysis in three stages. 
First, as discussed in Sj3l we modeled the continuum alone 
to determine the range of process parameters (Td, tr) con- 
sistent with the continuum light curve. We use this dis- 
tribution of estimated Td and a as a prior for the joint 
models of the continuum and line light curves in order to 
avoid the secondary solutions with Td — > as discussed 
in Sj2j Second, for each joint model, we find the best-fit 
top hat transfer function fEquation l22[) which maximizes 
the model likelihood calculated by Equation (fT?) . along 
with an updated set of process parameters. Finally, we 
ran an MCMC analysis on each joint model to calculate 
the statistical confidence limits on each best-fit param- 
eter found by global optimization on a grid, especially 
the time lag. We then compare these estimates to those 
derived from previous CCF analyses. We refer to these 
models as the "single- line" fits since they are solving for a 
single top-hat transfer function. The dotted histogram in 
Figure Q] shows the x 2 /dof distribution of the single- line 
model. It has a similar shape to the x 2 /dof distribu- 
tion of the stochastic model for only the continuum light 
curve, and confirms that the statistical model provides a 
good fit to the quasar variability, as well as the overall 
ansatz that the H/3 variability is a scaled and smoothed 
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Fig. 6. — Comparison between independent (single-line) and joint (two-line) fits to the Ho and H/3 light curves of NGC 3516. The 
red solid lines are the estimate from the single-line fits, while the blue dashed lines are those from the two-line fits. The top left (right) 
panel compares the likelihood distributions of the two fits for the Ha (H/3) line. The interval between the two dotted lines corresponds 
to a 3a range in the likelihood, while the two blocks above indicate the ±lcr range of the CCF peak analysis (upper) and CC centroid 
distribution (lower), where the central lines mark the T pea k and r cen t values, respectively. The two bottom left panels shows the color-coded 
covariance map between the two lags for the single-line and two-line fits, respectively. The contours in the bottom right panel compare 
the confidence levels calculated from MCMC for the Hcr/H/3 lags near the peak (black boxes inside left two panels). Working outward, the 
three contour curves are for la, 2a and 3a levels, respectively. Note that those are all observed-frame lags and the H/3 light curve here is 
the older of the two we have for NGC 3516. 



version of the continuum. The x 2 /dof distribution of 
the single-line model is somewhat worse than for fitting 
the continuum alone, but still reasonably consistent with 
statistical expectations. 

For the sake of uniformity of emission-line species in 
the comparison between the SPEAR and the CCF meth- 
ods, and to avoid confusion in the figures for sources with 
multiple line observations, we will focus on the the 66 
H/3 light curves in our subsequent analyses, and tabulate 
all the other emission-line lags we successfully computed 
with SPEAR in Table [U 

Figure [4] shows the comparison between CCF centroid 
time lags tccf and our lags tspear for all the H/3 lines. 
The range of uncertainties for tccf contains 68.3% of 
Monte Carlo realizations in the cross-correlation centroid 
distribution (CCCD), while our estimated error bound- 
aries are defined by the 68.3% (bxC/£ max = 0.5) con- 
fidence levels that encloses the best-fit lags (i.e., ±ler 
errors if the probability distribution is Gaussian in both 
cases) . Based on the structure of the lag probability dis- 
tribution, we can divide the "single-line" fits into five 
quality groups: 

(I) In most of the cases (43 of 66), the likelihood dis- 
tribution for the lags has a single peak and there is 
an unambiguous H/3 lag. 

(II) In 9 cases, the likelihood distribution has multi- 



ple peaks with significant (> 3er) likelihood differ- 
ences. This occurred for one season of Akn 120 
(JD49980-50175fl Mrk 110 (JD48953-49149), and 
Mrk 590 (JD49183-49338); two seasons of Mrk 79 
(JD48193-48393 and JD49996-50220); NGC 4051, 
PG 084^3, PG 1411, and PG 1617. Compared to 
our estimate, the CCF analysis picks a lower like- 
lihood peak or aliases several peaks into one broad 
peak. Generally, the two peaks are so close that the 
differences between the results from the two meth- 
ods are insignificant compared to the uncertainties. 

(Ill) In 7 cases, the likelihood distribution has multiple 
peaks of comparable significance (< 3cr: one series 
of NGC 3516 (JD47894-48047), Fairall 9, PG 0026, 
PG 0052, PG 1211, PG 1226, and PG 1307). They 
are shown in Figure 0] as the objects with a dashed 
line connecting the possible solutions. The tradi- 
tional CCF method seems to find one broad peak 
for these sources, rather than multiple peaks, lead- 
ing to large reported uncertainties for the estimate 
of tccf- These degeneracies are largely caused 

5 For brevity, we retain only the five least significant digits of 
the Julian Date. 

6 For brevity, we truncate the PG coordinate identifiers to right 
ascension only since this introduces no ambiguity in the present 
small sample. 
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Fig. 7. — Comparison between independent (single-line) and joint (two-line) fits to the Ha and H/3 light curves of PG 0026. The format 
is the same as in Figure [6] 



by poor light curve sampling that allows the light 
curve of the emission line to be mapped into the 
sampling gaps of the continuum. This problem is 
worst for the PG objects, which have many "sea- 
sonal gaps" over the long observing baselines (~ 
7.5 yr), leading to a clustering of solutions around 
180 days in the observed-frame. Such seasonal 
aliasing problems affect the CCF-based methods 
as well (jGrier et al.ll2008[) . 

(IV) In four cases, the light curves are very poorly sam- 
pled: IC 4329A, one season on NGC 4593 
(JD47894-48049), one season of Mrk 279 
(JD47205-47360), and one season of NGC 3227 
(JD47894-48 045). These cas e s were also flagged as 
unreliable bv lPeterson et abl <|2004f ) . so we exclude 
them from our subsequent analyses. 

(V) The lags derived from the SPEAR method appear 
to be wrong in two cases, 3C 120 and PG 1613. We 
also exclude both from our subsequent analyses. 
The 3C 120 light curves have a baseline of 7 years, 
but are very sparsely sampled. The CCF method 
finds a lag ~ 40 days in the observed frame. Al- 
though we find a sub-peak at 40 days, the model 
favors another peak of much higher significance at 
259 days. For PG 1613 we obtain a lag of ~ 575 
days in the observed-frame, much larger than the 
~ 50 day CCF estimate. In both cases, the longer 
lag is favored because it minimizes the data over- 
lap — 259 and 575 days put most of the line data 
in the seasonal gaps, and many points also lie be- 
fore the start of the continuum light curve. This 



is essentially an aliasing problem in our method. 
We also note that the continuum flux varied by up 
to 50% over the 7 year span of the light curves. 
We know empirically that the scaling coefficient A 
in the transfer function is inversely correlated with 
ionizing continuum flux (see the right panel of Fig- 
ure [10] and the discussion in $8]) , but we treat A as 
a constant parameter in each individual fit. This 
may create problems for light curves with the sig- 
nificant long term trends observed for these objects. 
Allowing A to vary and adopting a prior that pe- 
nalizes large lags that minimize light curve overlap 
would likely solve these problems. 

5. MODEL TEST FOR CORRELATED ERRORS 

Correlated errors have long been viewed as a prob- 
lem in traditional CCF analysis. Observations made at 
a common epoch are inevitably correlated by the pro- 
cesses required for calibration, light curve extraction, 
broad/narrow line modeling and removal of host or Fen 
contamination. Because no assumption about the prop- 
erties of the noise matrix N was made in Sj2j it is easy 
to include the effects of correlated errors within our ap- 
proach. While we did not make an extensive survey of 
our ability to model noise correlations between the con- 
tinuum and lines, we did carry out some experiments for 
objects noted a s pote ntially having strong covariances by 
iPeterson et"aT1 (|2004h . 

The simplest test is to introduce a covariance factor 
— 1 < r < 1 and add off-diagonal terms to the noise ma- 
trix ./V for line and continuum points measured at the 
same epoch of N c i(t,t) = Ni c (t,t) = rai(t)a c (t) in order 
to examine the sensitivity of the lag estimates to corre- 
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lated noise between the line and the continuum measured 
at each epoch. This should be present in the data at some 
level because of the challenge of consistently subtracting 
the contribution of the host galaxy to the line and the 
continuum in the presence of variable observational con- 
ditions. 

Figure [5] illustrates the effects of adding off-diagonal 
correlated noise terms on the H/3 lag estimate of PG 0844. 
The shift in the estimated lag (left top panel) induced by 
r varying from —1 to +1 is only about 0.25 days, much 
smaller than the median sampling interval of the light 
curves. The corresponding change in the likelihood (left 
bottom panel) shows a plateau at r > and slowly 
asymptotes to a maximum at r = +1, suggesting that the 
errors in the two light curves are positively correlated. 
The lag likelihood distribution (right panel) changes if 
we assume different levels of correlations r between the 
light curves. While the overall lag likelihood is greatly 
depressed in the r = — 1 case, the likelihood distributions 
are nearly identical in the r = and r = +1 cases. How- 
ever, the peaks near the best lag estimate (~ 12 days) 
arc slightly more significant in the r — +1 case than in 
the r = case. We explored this issue for several other 
systems, and generally the impact on the estimated lag is 
negligible, although different levels of (anti-)correlations 
are detected. 

6. JOINT ANALYSIS OF MULTIPLE LINES 

In 21 we found that poor light curve sampling was a 
significant problem in many systems, particularly in ob- 
jects with observed-frame lags on time scales similar to 
the seasonal gap spacing. However, if multiple lines have 
been measured, then we have significant, additional data 
to better sample the light curves under our overall ansatz 
that all light curves are scaled, smoothed, and displaced 
versions of the continuum. Simultaneous fits also de- 
termine the covariancc between the lags of the different 
lines. In this section, we explore simultaneously fitting 
the continuum and two emission-line light curves (here- 
after "two-line" fits, as opposed to the "single-line" fits in 
21 as we are now fitting two top- hat transfer functions) . 

Figure [6] summarizes the significant improvement in 
estimating the H/3 time lag of NGC 3516 of the 
IWanders et all (fl993h data (JD47894-48047) after in- 
cluding the Ha light curve (two-line) compared to us- 
ing the H/3 line alone (single-line). NGC 3516 is a case 
where the single-line H/3 fits shows a secondary peak at 
~ 42 days whose likelihood relative to the main peak 
at ~ 6 days is high, ln(£2nd / ' C-max) = —1-5 (solid curve 
in panel b). The Ha fit does not show such a secondary 
peak (panel a) . When we fit both simultaneously, the Ha 
light curve together with its well-determined lag adds ex- 
tra information to the continuum light curve, and thus 
better constrains the H/3 lag. The second H/3 peak is 
suppressed and there is a single unambiguous H/3 lag for 
the two-line fit (dotted curve in panel b). The improve- 
ment is most clearly seen in the structure of the Ha /H/3 
lag likelihood plane (panel c and d). If we zoom in on the 
remaining peak and run a MCMC chain using a flat prior 
on lags in the zoomed region, we can see that the two- 
line fits not only suppress the secondary peaks but also 
shrink the uncertainties in the primary peak to produce 
better results for both lines (panel e). 

The joint analysis of multiple lines is especially useful 



for the PG objects, whose light curves show observational 
gaps of period ^180 days in the observed-frame. In the 
single line fits, the model would always show (sub)pcaks 
for lags ^180 days because of the seasonal aliases (the 
seasonal stripes in Figure HJ). It is not possible, however, 
to do this for 2 lines simultaneously, so the two-line fits 
largely eliminate seasonal aliasing. Figure [7] illustrates 
this for PG 0026. In particular, the broad H/3 likelihood 
distribution shrinks significantly and the maximum like- 
lihood lag drops from ^160 days to ^106 days (~ 140 to 
~ 93 in the rest frame) and is in better agreement with 
the Ha results. Although the traditional CCF method 
makes similar estimates (green and blue bands in two 
top panels) , it yields significantly larger uncertainties by 
aliasing several peaks into one broad CCF centroid dis- 
tribution. 

We performed similar joint analyses for the 21 sources 
for which we have multiple emission line light curves and 
recompile the results for the H/3 lags, as shown in Fig- 
ure [8] Fortunately, all the sources whose H/3 lags were 
found to be ambiguous in the single- line fits (i.e., the 
7 H/3 lags from groups III in 2} are improved by the 
two-line fits, although the degree of improvement varies. 
We also dropped lag estimates that were either flagged 
as unreliable or believed to be wrong (i.e., the 6 H/3 lags 
from groups IV and V in 21 an d keep only those objects 
deemed to give robust estimates of lag by our method 
(i.e., the 60 H/3 lags from groups I, II and III). To illus- 
trate the quality of the final result for each source, we 
divide all 60 remaining sources into 4 new groups based 
on the results of both the single-line fits in § 0] and the 
two-line fits, using different symbols for the 4 new groups 
in Figure El 

(A) The 43 group I light curves from § @] with a single 
unambiguous H/3 lag. Seven of the objects have 
light curves of lines other than H/3 to carry out 
two-line fits, but they provided little gain when the 
single- line fits already provided good lag estimates. 

(B) The 10 group II sources from § [4] with a robust 
H/3 lag estimate but potentially larger uncertainties 
due to the presence of low significance (> 3tr) sub- 
peaks in the lag likelihood distribution. Most of 
those sources do not have the multiple line light 
curves needed to carry out two-line fits. 

(C) The four group III sources (NGC 3516, PG 0026, 
PG 0052, and PG 1307) from g] with multiple 
peaks in the single-line lag likelihood distribution 
where the ambiguity is removed by the two- line fits. 

(D) The three group III sources (Fairall 9, PG 1211, 
and PG 1226) from 21 with multiple peaks in 
the single-line lag likelihood distribution where the 
two-line fits fail to remove the ambiguity. We 
picked the most significant peak as the solution and 
extended the uncertainty to cover all the possible 
solutions. 

Recall that we have dropped the 6 group IV and V 
sources (IC 4329A, NGC 4593, one season of Mrk 279, 
one season of NGC 3227, 3C 120, and PG 1613) out 
of all 66 H/3 light curves following the discussion in 21 
Green circles, blue pentagons, dark violet squares, and 
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Fig. 8. — Comparison of the H/3 time lags from CCF analysis and the SPEAR method, similar to Figure[4] updated where we have used 
the two-line fits and dropping the 6 unreliable sources. Four types of symbols are used to indicate our estimate for increasing levels of 
ambiguity in the lag estimate. Objects with inconsistent lag estimates between the two methods are labeled. 



red triangles correspond to sources of group A, B, C, and 
D, respectively. There is general agreement between the 
two methods, but also several discrepancies, as 7 of our 
H/3 lag estimates are inconsistent with the CCF results 
given their error estimates. We marked these sources in 
Figure [5] and now discuss each case individually, 

NGC 7469.— We estimate an H/3 lag of ll-T^ days, 

as opposed to tccf =4.7l^s- However, if we use a 
Dirac delta function for the transfer function instead of 
a tophat, the estimated time lag changes to 4.3 days, 
in agreement with the CCF result. Thus, the discrep- 
ancy originates from the improvement of fit with a tophat 
smoothing kernel. The continuum of NGC 7469 was in- 
tensively monitored to searc h for time lags betw een the 
UV and optical continuum (|Collier et al.l 119981 ). so its 
continuum light curve is densely sampled while the H/3 
light curve is much less so. The model has to smooth 
the continuum light curve heavily (i.e., a broad tophat 



width) to obtain a good fit, which at the same time shifts 
the time lag estimate to a longer value than it would be 
with a zero width (i.e., a delta function). This is sugges- 
tive of the continuum errors being underestimated, or a 
more realistic transfer function is required. 

Mrk 79 (years 2 and 4). — In both cases, we estimate 
larger time lags than the CCF results, although there are 
sub-peaks which correspond to the CCF lags. For year 
2 (JD47838-48044), while the CCF centroid gives a lag 
of 16.4lg 7 days, the CCF peak estimate is 19**2 days, 
more consistent with our estimate of 30. 91 1 , \ days. For 
year 4 (JD49996-50220). iPeterson et al.l (l200l flagged 
it as "unreliable" for the poor light curve sampling. Our 
method shows a dense array of sub-peaks in the lag like- 
lihood distribution, but the most significant peak is at 
43.6+^ days. 
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Fig. 9. — The Rblr-L relation for H/3. The luminosity is AL>,(5100A) and the BLR radius is equivalent to the lag in units of light 
days. The open symbols and gray solid circles indicate the measurement from SPEAR method and from CCF method for the same set of 
sources, respectively. The gray solid curve is the fit to the CCF Bblr~^ relation, while the rest of the curves are the fits to the SPEAR 
-Rblr - ^ relation, using four subsets of the sources (see Table [2] for details of each fit). The slope of the fit to the SPEAR Rblr~L relation 
a is steeper than the CCF relation, but the two are consistent within the uncertainties. cr rms is the rms scatter of each fit. 



PG 0844. — As discussed in § the CCF estimate of 
the H/3 lag (34.4i^'j) for PG 0844 is likely susceptible 
to correlated errors, while our method estimates a lag of 



lag, favoring a longer lag that is more consistent with 
lags of the other Balmer lines. 



12.2_ 1 ; 3 days regardless of the value of correlation coef- 
ficient r. 



We estimate an H/3 lag of 149.3l*'g days, 



PG 1426.— We estimate an H/3 lag of 161.6 
as opposed to tccf = 103.2 



+32.5 
-40.3- 



±ii 9 i days, 
Similar to PG 0052 



PG 0052. 

as opposed to tccf ^^-27 t- The single-line fit shows 
multiple peaks and usually one would be inclined to mis- 
trust a peak at the seasonal alias (a rest-frame lag of 
150 days corresponds to 170 days in the observed- frame). 
However, the joint Ha/H/3 fit clearly reinforced this so- 
lution. 



PG 1307.— We estimate an H/3 lag of 188.8 
as opposed to t C cf = 121.9is3;|. The j° int Ha/H/3 fit 



+ 3 7 7 days, 



suppressed the false peak which corresponds to the tccf 



and PG 1307, the joint Ha/H/3 fit reinforced a solution 
which is otherwise susceptible to the seasonal gap effect. 

We carried out a similar analysis for each data set, in- 
cluding all emission lines besides H/3, as summarized in 
Table [T] Note that in the table we only include 87 light 
curves for which we have successfully computed lags. The 
object is identified in column (1). The emission line 
and its light curve Heliocentric Julian Date range are 
listed in columns (2) and (3), respectively. Column (4) 
gives the rest-frame time lag estimate from the SPEAR 
method, while column (5) indicates the associated "am- 
biguity" (i.e., the group membership) defined above. 
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Fig. 10. — Lag (left) and scaling coefficient (right) of the H/3 transfer function as a function of continuum luminosity from 14 years of 
NGC 5548 data. The best-fit slopes are also reported for each panel and shown by the black solid lines. The black dashed line in the left 
panel is the best fit with a fixed slope of 0.50. The number inside each solid circle indicates the year of observation for each light curve 
starting from Dec. 1988. Note that in the left panel the point for year 12 is hidden under that of year 13. 



7. THE Bblr-L RELATION FOR H/3 

With the revised set of H/3 lags, and the starlight- 
corrected optical luminosity of each AGN from 
IBentz et all (pOOl) . we have calculated the fit to the 
Rblr~L relationship for our sample 

i?BLR = C '" ( 1043-5 CTgss-l) (23) 

and compared it to that based on CCF lags in Figure [9l 
We obtained a slope a = 0.579±0.010 for all the SPEAR 
lags regardless of the level of "ambiguity" at which we 
probe the slope. This slope is slightly steeper than pre- 
vious estimates, and only marginally consistent with the 
naive theoretical prediction of a = 0.5. Compared to the 
CCF-based -Rblr - -t' fit of the same sample of AGNs (blue 
filled circles), our Rblr~L fit has a steeper slope but 
comparable rms scatter cr rms , which grows smaller as we 
use more reliable lags. Table [2] gives the results from the 
different fits using the 4 combinations of groups indicated 
by column (1). Column (2) gives the number of data 
points used in each fit. We fit each combinatorial data 
set using lag estimates from both the SPEAR (columns 
3-6) and CCF methods (columns 7-10). Two param- 
eters in Equation (|23[) are listed in columns (3) and 
(4) for SPEAR method, and in column (7) and (8) for 
CCF method, respectively. Column (5) and (6) give the 
X 2 /dof and rms scatter for our fit, while column (9) and 
(10) give these statistics correspondingly for the CCF 
method. 

Our Rblr-L fits have a larger \ 2 /dof than the CCF 
ones. This does not necessarily mean they are poorer fits, 
because our lag estimates generally have tighter error- 
bars than the CCF estimates. It could indicate that our 
approach underestimates uncertainties, that the CCF 
method overestimates uncertainties, or that we are not 
taking into account intrinsic scatter in the Rblr~L re- 
lationship. Since most of the group C and D sources 



are high-redshift luminous PG objects, the rms scatter 
for our method decreases from 0.229 dex to 0.196 dex 
after dropping them from the fit. Three outliers from 
the CCF Rblr-L relation (NGC 7469, years 1 and 4 
of Mrk 79) are also the sources where our lag estimates 
are inconsistent with CCF results. When we use our lag 
estimates, these three CCF outliers lie on the Rblr-L re- 
lation, which reduces the rms scatter near L opt ~ 10 43 5 
ergs s _1 . Note that there is significant scatter in the 
Rblr~L relation even for multiple estimates for a sin- 
gle source, as shown in the left panel of Figure [TU] for 
NGC 5548. 

8. Rblk-L RELATION OF NGC 5548 REVISITED 

So far, we have carried out our calculations assuming 
that the parameters are constant during a season. This 
is likely true for the underlying variability process. If we 
model either the full continuum light curve or the individ- 
ual seasons, we find estimates for the process parameters 
Td and a that are statistically consistent. We do observe 
lags that vary from season to season, and these are ar- 
guably correlated with luminosity. If so, they should also 
be varying within seasons, and we have not accounted 
for this. Similarly, we assume the scaling between the 
continuum and line fluxes does not vary over a season, 
although we do observe it to vary between seasons. 

The nearby Seyfert 1 galaxy NGC 5548, with its many 
continuous years of monitoring data, serves as an ideal 
example of an AGN changing its variability levels from 
season to season. Figure [TUJ illustrates the continuum 
flux dependence of both the H/3 lag (r) and the scal- 
ing coefficient A for 14 seasons of NGC 5548 data. We 
clearly see trends that the lag increases with luminosity 
and the amplitude of the response diminishes. If we fit 
the lag, we find a steep slope, (r) cx £°- 73 ± - 10 that is 
inconsistent with the expected (t) oc L 5 . However, the 
poor fit (x 2 /dof = 4.17) suggests that either the uncer- 
tainties arc underestimated or intrinsic scatter dominates 
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TABLE 1 
Rest-Frame Lag Estimates 



Julian Dates ^spear 



Object 


Line 


(-2400000) 


(days) 


Group 


(1) 


(2) 


(3) 


(4) 


(5) 


3C 390.3 


H/3 


49718- 


-50012 


27.91?-* 


A 


3C 390.3 


Lye? 


49718- 


-50147 


11.9lf 6 5 


D 


3C 390.3 


Civ A1549 


49718- 


-50147 


15.0l 2 3 °o 


C 


Akn 120 


H/3 


48148- 


-48344 


35.7l 6 9 ; 7 2 


A 


Akn 120 


H/3 


49980- 


-50175 


29.7+fl 


B 


Fairall 9 


H/3 


50473- 


-50665 


i9.4±*y 


D 


Fairall 9 


Hen A1640 


50473- 


-50713 


12.0+™ 


C 


Fairall 9 


Lya 


50473- 


-50713 


in 1+0.5 
— 0.5 


C 


Mrk 79 


H/3 


47838- 


-48044 


25. 5_ 14 4 


B 


Mrk 79 


H/3 


48193- 


-48393 


30.91J* 


A 


Mrk 79 


H/3 


48905- 


-49135 


n.2+ 7 2 % 


B 


Mrk 79 


H/3 


49996- 


-50220 


43.6l^ 


A 


Mrk 110 


H/3 


48953- 


-49149 


rjr + 2.3 

13 1 


B 


Mrk 110 


H/3 


49751- 


-49874 


33.9l 6 5 ;'i 


A 


Mrk 110 


H/3 


50010- 


-50262 


21. 51 2 ; 2 


A 


Mrk 279 


H/3 


50095- 


-50289 


18. 3+5'j 


A 


Mrk 290 


H/3 


54184- 


-54301 


- -+0.7 

■ —0 5 


A 


Mrk 335 


H/3 


49156- 


-49338 


15.3111 


A 


Mrk 335 


H/3 


49889- 


-50118 


12.9ll;« 


A 


Mrk 509 


H/3 


47653- 


-50374 


69.9 + °| 


A 


Mrk 509 


Hen A4686 


47653- 


-50374 


52.2l^ 


D 


Mrk 590 


H/3 


48090- 


-48323 


19.0l 2 j! 


A 


Mrk 590 


H/3 


48848- 


-49048 


19.5l 2 4 .° 


A 


Mrk 590 


H/3 


49183- 


-49338 


32.6l|| 


B 


Mrk 590 


H/3 


49958- 


-50122 


30.9l2^ 


A 


Mrk 817 


H/3 


49000- 


-49212 


20.9l 2 | 


A 


Mrk 817 


H/3 


49404- 


-49528 


17.212,? 


A 


Mrk 817 


H/3 


49752- 


-49924 


35.9l*; 8 8 


A 


Mrk 817 


H/3 


54200- 


-54330 


10.8l^ 


A 


NGC 3227 


H/3 


48623- 


-48776 


10.6l^ 


A 


NGC 3227 


H/3 


54180- 


-54273 


4.4l° - 3 5 


A 


NGC 3516 


Ha 


47894- 


-48047 


14.0l° - 7 7 


A 


NGC 3516 


H/3 


47894- 


-48047 


6-ll°o' 5 7 


C 


NGC 3516 


H/3 


54181- 


-54300 


14.6111 


A 


NGC 3783 


H/3 


48607- 


-48833 


? O + 0.3 

' -°_o 7 


A 


NGC 4051 


H/3 


54180- 


-54311 


2.5l° i 


B 


NGC 4151 


H/3 


53430- 


-53471 


6 +0 - 6 


A 


NGC 4593 


H/3 


53430- 


-53471 


4 5 +0 -' 7 


A 


NGC 7469 


H/3 


50237- 


-50295 


11.7l° - 5 7 


A 


NGC 7469 


Si iv A1400 


50245- 


-50293 


2.0l° * 


A 


NGC 7469 


CivA1549 


50245- 


-50293 


10.6 + S" 2 2 


A 


NGC 7469 


Hen A1640 


50245- 


-50293 


o.8l° ; 2 2 


A 


PG 0026+129 


Ha 


48836- 


-51084 


88-Olsl 


B 


PG 0026+129 


H/3 


48545- 


-51084 


92.7l 7 ;° 6 


C 


PG 0052+251 


Ha 


48837- 


-51084 


157.6l 2 j 


A 


PG 0052+251 


H/3 


48461- 


-51084 


149.3+ 42 


C 


PG 0052+251 




48461- 


-51084 


154.911J 


C 


PG 0804+761 


Ha 


48319- 


-51085 


133.4l|J 


C 


PG 0804+761 


H/3 


48319- 


-51085 


116.8l 2 3 


A 


PG 0804+761 


H-j 


48319- 


-51085 


71.1+* 6 i 6 


D 


PG 0844+349 


Ha 


48319- 


-51085 


20.8l°* 


A 


PG 0844+349 


H/3 


48319- 


-51085 


12.2+j- 2 


B 


PG 0844+349 




48319- 


-51085 


17 7 +2.-5 
— 2.1 


C 


PG 0953+414 


H/3 


48319- 


-50997 


162.2l 2 j 


A 


PG 0953+414 


H-; 


48319- 


-50997 


1 60 2+ 3 3 


D 


PG 1211 + 143 


Ha 


48319- 


-51000 


76.3l° ; 7 5 


C 


PG 1211+143 


H/3 


48319- 


-51000 


7 „ 0+0.9 

loo -25.4 


D 


PG 1211 + 143 


H-; 


48319- 


-51000 


"'■'-10.1 


A 


PG 1226+023 


Ha 


48361- 


-50997 


380.0 + f 7 


B 


PG 1226+023 


H/3 


48361- 


-50997 


105.5+ 28 2 4 1 


D 


PG 1226+023 




48361- 


-50997 


263.8l« ° 6 


A 


PG 1229+204 


Ha 


48319- 


-50997 


45-7l?;i 


A 


PG 1229+204 


H/3 


48319- 


-50997 


42.8l 2 i 


A 



TABLE 1 
— Continued 



Object 

(i) 


Line 
(2) 


Julian Dates 
(-2400000) 
(3) 


1-SPEAR 

(days) 
(4) 


Group 

(5) 


PG 1307+085 


Ho 


49130-51000 


189 


A 


PG 1307+085 


H/3 


48319-51042 


188. 8 + 3 7 


c 


PG 1307+085 


H-T 


48319-51042 


218 9 + I; 2 , „ 

— 124.8 


D 


PG 1411+442 


Hn 


48319-51038 


59 3+1 , 1 


A 


PG 1411+442 


H/3 


48319-51038 


53 5+i 3 , 1 


B 


PG 1426+015 


H/3 


48334-51042 


161 6 + ?; 9 , 
u u — 11.1 


A 


PG 1617+175 


Ha 


48362-51085 


106.9 + ^3 S 3 


B 


PG 1617+175 


H/3 


48362-51085 


88 2+ 31 „° 

uu -5.9 


B 


PG 2130+099 


H/3 


54352-54450 


23 2 + t'i 


A 


PG 2130+099 


He 11 A4686 


54352-54450 


32 + 

w — 4.5 


A 


NGC 5548 


H/3 


47509-47809 


21 2+°-» 
— 1.0 


A 


NGC 5548 


H/3 


47861-48179 


W.3 + °l 


A 


NGC 5548 


H/3 


48225-48534 


15 8 +2 -'; 


A 


NGC 5548 


H/3 


48623-48898 


"•oil: 2 , 


A 


NGC 5548 


H/3 


48954-49255 


J.O.O_3 Q 


A 


NGC 5548 


H/3 


49309-49636 


10.811* 


A 


NGC 5548 


H/3 


49679-50008 


24.2i£;i 


A 


NGC 5548 


H/3 


50044-50373 


i6.il° ; 3 6 


A 


NGC 5548 


H/3 


50434-50729 


l6.8l° ;* 2 


A 


NGC 5548 


H/3 


50775-51085 


26.9l 2 ; 2 


A 


NGC 5548 


H/3 


51142-51456 


23.8l 2 j 


A 


NGC 5548 


H/3 


51517-51791 


H H + 1 - 3 
S °-3.9 


A 


NGC 5548 


H/3 


51878-52174 


8.7l° ; 5 5 


B 


NGC 5548 


H/3 


54180-54332 


l6.3+i;» 


A 


Note. — Lag estimates and 


confidence limits for Groups 


A and B 



are calculated by the single-line fits, while those for Groups C and 
D are from the two-line fits. 

the goodness-of-fit. If we rescale the uncertainties so that 
the best-fit model has % 2 /do/ = 1, the flatter L 5 slope 
is not ruled out, with a A\ 2 of only 0.16. 

These problems can be addressed by making the lags 
and the line-to-continuum scaling a function of the con- 
tinuum luminosity. For the luminosity dependence of 
lags, the simplest approach would be to de-lag the line 
light curve as (r) cx L a instead of shifting the entire light 
curve by the same ti ag , and then optimize the fits over 
the additional parameter a. Unfortunately, we cannot 
fit the full NGC 5548 light curve because the resulting 
matrix dimensions are impractically high (_ftT=3085 data 
points). We instead estimate the normalized likelihood 
distribution for a in each season and then combined the 
likelihoods, as shown in Figure [TT] (we did not include 
year 13, which was part of the less reliable group B). 
This "breathing" effect is clearly detected, and the log- 
arithmic slope estimate of a = 0.44lQ Qg is consistent 
with the naive expectation a = 1/2 and the Rblr~L re- 
lation in Figure HI Using almost the same set of light 
curves fr om NGC 5548 (we add year 14 and exclude 
year 13). lCackett fc Hornd (|2006h find a much shallower 
slope (0.1-0.46) with a luminosity-dependent delay map, 
in better agreement with the prediction of p hotoioniza- 
tion models (~ 0.23; iKorista fc Goadll2004f) . However, 
their small correction for the host galaxy s tarlight may 
artifi cially flatten their estimate of the slope (|Bentz et al.l 
l2007t ). Note that for this experiment we did not make 
the linc-to-continuum scaling coefficient A a function of 
continuum luminosity in the fit. Such a full scale calcu- 
lation should be carried out using the complete data set. 
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TABLE 2 
BLR Size-Luminosity Relation 



Groups 


N 


CSPEAR 


<*SPEAR 


^SPEAR 


SPEAR 
rms 


CCCF 


"CCF 


XCCF 


CT CCF 
rms 


Included 




(lt-days) 




(dex) 


(lt-days) 




(dcx) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


A,B,C,D 


60 


1.36 ±0.01 


0.579 ±0.010 


8.13 


0.229 


1.32 ± 0.01 


0.513 ±0.017 


3.50 


0.203 


A,B,C 


57 


1.36 ±0.01 


0.580 ±0.011 


8.52 


0.207 


1.33 ±0.01 


0.519 ±0.019 


3.54 


0.205 


A,B 


53 


1.36 ± 0.01 


0.558 ±0.013 


8.39 


0.206 


1.33 ± 0.01 


0.521 ±0.020 


3.81 


0.213 


A 


43 


1.36 ±0.01 


0.556 ±0.013 


9.44 


0.196 


1.32 ±0.01 


0.518 ± 0.020 


4.22 


0.211 



a= 0.48ig|° (±lcr CL) NGC 5548 




0.2 0.4 0.6 0.8 

a 



Fig. 11. — Likelihood distribution of a for 13 years of NGC 5548 
light curves. The normalized log-likelihood is calculated by adding 
the likelihood distribution functions for the 13 individual years 
together. 

9. DISCUSSION 

We have demonstrated that direct fitting of contin- 
uum and line light curves is a viable approach to mea- 
suring reverberation lags , confirming the initial study of 
iRybicki fc Kleviial ()1994l ). It provides a full statistical 
framework for determining time lags and estimating their 
uncertainties, including the full contributions from cor- 
related noise, de-trending and interpolation. In essence, 
the lags are determined using a weighted average of all 
statistically acceptable models for interpolating the un- 
derlying true light curve. While we used the assump- 
tion that the underlying variable process had an expo- 
nential correlation function corresponding to a damped 
random walk, any other statistical p rocess could be sub - 
stituted. We note, h owev er, that [ Kelly et ail (120091 ) . 
iKozlowski et all (|2010t) and iMacLeod et al.l (|201dh have 
found the exponential correlation function to be an excel- 
lent model of quasar light curves, just as we have found 
here, although we modeled the light curves in flux rather 
than magnitude. 

Because we are explicitly modeling the light curves, we 
must include an explicit model of the transfer function. 
Here we used a top hat for simplicity. It includes the sim- 
ple limits of a delta function and a uniform thin shell, and 
is likely a reasonable model for any sin gle-peaked transfer 
funct ion given the available data (see lRvbicki fc Klevnal 
I1994D . As with the model for the variability process, us- 
ing an alternative transfer function simply requires com- 
puting the appropriate terms of the covariance matrix. 



Aside from the case of NGC 7469 where it seemed to 
affect the lag estimation, we did not discuss the tophat 
width. In general, there is a relatively strong degeneracy 
between At, the width, and A, the scaling between the 
continuum and line light curves. When At is large and 
the continuum is heavily smoothed, the model will try to 
increase the variability amplitude by artificially boosting 
A to re-align the continuum and line light curves. How- 
ever, the degeneracy does not seem to lead to problems 
in estimating the mean lag unless the line light curve is 
very poorly sampled. The traditional CCF method does 
not implicitly assume a shape for the transfer function 
but calculates the lag as either the barycenter (r cent ) 
or the peak (r pea fe) of underlying transfer function (con- 
volved with data). The difference between the two some- 
times can be large and hard to reconcile unless the trans- 
fer function can be modeled explicitly. For future high- 
fidelity datasets, our approach should also have no diffi- 
culty constraining the shape of transfer functions. 

The most important future path for this method is to 
simultaneously fit multiple line components, whether dif- 
ferent lines (e.g., H/3, Ha, etc.), velocity sub-components 
of individual lines or multiple continuum bands. As long 
as the overall ansatz that all light curves are scaled and 
smoothed versions of the continuum holds, combining 
many light curves with differing lags means that the lag 
estimate for any given light curve is now derived from 
a better sampled estimate of the continuum variabil- 
ity. A second advantage, particularly for attempts to 
study the velocity structure of a particular broad line, 
is that such joint analyses will correctly infer the co- 
variances between the individual lags. Current velocity- 
dependent lags have uncerta i nties comparable to their 
differ ences (|Bentz et al.ll200l l2010l : iDennev et aDl200l 
I20101 ). but it may be true that these differences actu- 
ally have a strong covariances, so that the differences 
are far more significant than estimates from analyzing 
the light curves in isolation. The method can also allow 
for luminosity-dependent lags or line-continuum scaling 
factors. Also note that while we only use the linear pa- 
rameters of the model to remove the light curve means, it 
is a very flexible tool for de-trending or cross-calibrating 
light curves whose model uncertainties will be fully in- 
cluded in lag estimate. 

The most important observational implication of this 
approach is the value of measuring multiple lines, es- 
pecially those with high ionization potentials. In our 
approach, multiple lines with differing lags allow one 
to overcome many of the sampling problems inherent 
to cross-correlation methodologies. At its simplest, one 
light curve can be aliased into a (seasonal) sampling gap, 
but two cannot be unless the transfer functions are sim- 
ilar (i.e., the lines have similar lags). Given the radial 



Measuring Reverberation Lags 



17 



ionization stratification of the BLR (Peterson 1991), the 
lag difference between two lines is proportional to the 
difference in their ionization levels. In this paper, how- 
ever, the lines we used for two-line fits are mostly pairs 
of two Balmer lines, which have similarly low ionization 
levels. Thus, the observational goal should be to obtain 
data for multiple lines with a broad range of ionization 
potentials. Indeed, with a wide variety of emission-line 
lines, it is in principle possible to combine the reverber- 
ation results with photoionization equilibrium modeling 
to highly constrain the geometr y and physics of the BLR 
(jHorne. Korista. k GoaJl200l . 

The only significant algorithmic challenge comes from 
the 0(K 3 ) scaling of the computational cost with the 
number of data points K. Unfortunately, the reverber- 
ation mapping problem is very different from simply us- 
ing the damped random walk to model the continuum 
light curves, where we can take advantage of the partic- 



ular structure of the covariance matrix to calculate the 
necessary matrix inversions in O(K) operations. Since 
the expensive matrix inversion is required for each like- 
lihood calculation, it becomes difficult to analyze large 
data sets, particularly if the number of parameters also 
increases greatly as in a full simultaneous model of lags 
as a function of line velocity. These problems can be ad- 
dressed using hyper-threaded or parallel versions of the 
underlying algorithm. 
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APPENDIX 

COVARIANCE MATRIX OF THE CORRELATION FUNCTIONS 

The expressions for the covariance matrices used in this paper and the accompanying code assume that the transfer 
function is a simple top hat, 

= A(t 2 -t 1 y 1 for t 1 <t-t'<t 2 . (Al) 



For this transfer function, we can analytically calculate the correlation functions in Equation ©, © and 
tively. 



respec- 



1S 



The Covariance Matrix Between the Continuum and One Line 
The covariance between continuum s c (t) at tj and line si(t) at ti with transfer function defined as in Equation (|A1|) 

( e -t L /r d _ e -t H /r d if t L > o 

(*c(tj)*j(ti)) = T A CJ 2 A I e *«/Td _ e t L /r d if tH < o ( A2 ) 

[ 2 - e tL / Td - e~ tH ' TA iit L <0<t H , 

where ti, = ti — tj — t 2 and tn = U — tj — t\. 



The Covariance Matrix Between Two Lines 

Consider the case when the first line s/(i) has transfer function ^>(t — t 1 ) as defined in Equation (|A1[) and the other 
line s'i(t) has transfer function ^'(t — t') 

q>'{t-t') = BiU-tsy 1 for t 3 <t-t' <U, (A3) 
where t^ — t 3 <t 2 —t\. The covariance between line si (t) at time i, and line sj (t) at time tj (Equation [8]) is 



{s l (t i )s' l (t J ))=r2a 2 AB 



,-\tL\/r d 



\th\/r d 



where 



■h |/r d _ e 


-\tMl\/r d _ 


g— |*M2 


\/ T d ^_ 


( 2t H /r d iit M 2<0<t H 


{ 2(i 4 -i 3 )/rd i£t M 2<0<t H 










{-2t L /T d ift L <0<t M1 


H |/r d _ g 


-|*JUl|/Td _ 


g— |*M2 


l/Td 


if t L > or t H < 0, 


t L 


= (ti-tj) 


-(t 2 - 


h), 




tMl 


= (ti-tj) 


-(i 2 - 


t 3 ), 




t>U1 


= (ti-tj) 


-(h- 


t 3 ), 




and tu 


= (ti-tj) 


-(h- 


U). 





(A4) 



(A5) 

By definition, the covariance for the autocorrelation of line si(t) between time ti and tj (Equation[7]) can be obtained 
by equating *'(i - t') with *(t - t 1 ) so that B = A, t 3 = t x and t 4 = t 2 . 
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